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ABSTRACT 

The  ratio  estimation  problem  arises  in  many  different 
applications  settings.  This  paper  is  concerned  with  the 
interplay  between  gradient  estimation  and  ratio  estima¬ 
tion.  Given  unbiased  estimators  for  the  numerator  and 
the  denominator  of  the  ratio,  as  well  as  their  gradients, 
joint  central-l<mit  theorems  for  the  ratio  and  its  gradient 
are  derived.  The  resulting  confidence  regions  are  of  poten¬ 
tial  interest  when  optimizing  such  ratios  numerically,  or 
for  sensitivity  analysis  with  respect  to  parameters  whose 
exact  value  is  unknown.  The  paper  also  briefly  disensses 
low-bias  estimation  for  the  gradient  of  a  ratio. 

1.  INTRODUCTION 

Let  (A,  B)  be  a  pair  of  jointly  distributed  real-valued 
random  variables.  The  estimation  of  the  ratio  a  = 
E[A]/E[B ]  is  known,  in  the  simulation  literature,  as  the 
ratio  estimation  problem.  Such  ratio  estimation  problems 
arise  in  many  different  applications  settings.  For  example, 
it  is  well  known  that  the  steady-state  mean  of  a  positive 
recurrent  regenerative  stochastic  process  can  be  expressed 
as  such  a  ratio  of  expectations;  see,  for  example,  Section 
3.3.2  of  Bratley,  Fox,  and  Schrage  (1987)  or  Chapter  2  of 
WoUT  (1989).  In  Section  2  of  this  paper,  we  will  discuss  the 
ratio  estimation  problem  in  greater  detail  and  offer  addi¬ 
tional  examples.  It  will  turn  out  that  the  infinite-horizon 
discounted  coot  of  a  non-delayed  regenerative  process  can 
also  be  expressed  in  terms  of  an  appropriately  chosen  ra¬ 
tio  estimation  problenfc  This  fact  was  first  pointed  out  by 
Fox  and  Glynn  (1989). 

Recently,  the  simulation  community  has  devoted  a 
great  deal  of  attention  to  the  nse  of  simalation  as  an  opti¬ 
mization  tool.  An  important  component  of  this  research 
effort  has  been  the  development  of  estimation  method¬ 
ology  for  computing  the  gradient  of  a  real-valued  per¬ 
formance  measure  with  respect  to  a  (finite-dimensional) 
decision  parameter  vector.  Such  gradients  play  an  im¬ 
portant  role  in  many  iterative  algorithms  for  performing 
both  constrained  and  unconstrained  mathematical  opti¬ 
misation.  This  paper  is  intended  as  a  study  of  the  ques¬ 
tion  of  how  to  use  this  gradient  estimation  methodology 


in  the  setting  of  the  ratio  estimation  problem. 

The  paper  is  organized  as  follows.  In  Section  2,  a 
number  of  different  applications  in  which  ratio  estimation 
problems  arise  are  discussed,  and  the  mathematical  frame¬ 
work  for  the  remainder  of  the  paper  is  described.  Section 
3  is  devoted  to  deriving  a  confidence  interval  methodol¬ 
ogy  for  estimating  the  partial  derivative  of  a  ratio.  In  ad¬ 
dition,  a  joint  central-limit  theorem  for  the  simultaneous 
estimation  of  the  entire  gradient  is  obtained.  In  Section  4, 
low-bias  estimation  issues  are  discussed.  Finally,  Section  5 
discusses  some  experimental  results  related  to  gradient  es¬ 
timators  for  ratios,  and  Section  6  concludes  the  paper  with 
a  brief  summary.  The  proof  of  our  main  theorem  (Theo¬ 
rem  1)  is  given  in  the  Appendix.  The  other  proofs  are  not 
given  here.  A  (future)  more  elaborate  version  of  the  paper 
will  contain  all  the  proofs,  derive  a  joint  central-limit  the¬ 
orem  that  can  be  used  to  simultaneously  estimate  the  gra¬ 
dient  and  the  Hessian  of  mixed  second-partial  derivatives 
of  a  ratio,  and  provide  further  numerical  illustrations. 

2.  EXAMPLES  OF  RATIO  ESTIMATION  PROB¬ 
LEMS 

As  discussed  in  the  introduction,  the  ratio  estimation 
problem  is  concerned  with  the  estimation  of  the  ratio 


where  (A,  B)  is  a  pair  of  jointly  distributed  real-valued 
random  variables.  We  now  proceed  to  offer  several  exam¬ 
ples  of  this  estimation  problem. 

EXAMPLE  1.  Let  X  *  (X(t),  I  >  0}  be  a  real-valued 
(possibly)  delayed  regenerative  process  with  regenerative 
times  0  <  T(0)  <  T(l)  <  •  •  •.  Foe  t  >  1,  let 


Bi  m  T(0-T(i-1). 
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If  E[Ai  +  Bi]  <  oo,  then  it  can  be  shown  [sec,  for  example, 

Asmussen  1987,  or  Wolff  1989)  that 


lira  -  /  X(a)ds  =  a  =  £[A, ]/£[/?,]. 

1  Jo 

Hence,  as  discussed  in  the  introduction,  the  steady-state 
mean  of  such  a  process  can  be  expressed  as  the  ratio  of 
the  two  expectations  £[Ai]  and  E[B\\. 

EXAMPLE  2.  Let  X  =  (-V(t),  t  >  0}  be  a  non-delaved 
regenerative  process,  taking  \  abios  in  a  slate  space  5,  wit h 
regenerative  times  0  =  T(0)  <  T(l)  <  •  •  •.  Let  /  and  g 
be  two  real-valued  non-negative  (measurable)  functions 
defined  on  5,  and  set 

V{t)  =  /  g(X(a))da 

Jo 

a  =  exp{- V(t)]/(A'(t))dt  . 

Then,  a  is  the  infinite-horizon  expected  discounted  cost, 
the  process  g(X(t))  corresponds  to  the  (state-dependent) 
discount  rate  at  time  t,  and  /(A'(t))  is  the  (undisconnted) 
rate  at  which  cost  is  incurred  at  time  t.  A  common  choice 
for  g  is  the  one  in  which  </(•)  is  constant  and  equal  to 
p  >  0,  in  which  case 

o  =  E  |^°°  exP[— pt]/(A’(t))dt] 

is  the  infinite-horizon  p-discounted  cost.  Let 


A, 

Cy 

By 


exp 


g(X(s))<h 


exp  [-VTT(l))] 


f(X(t))dt 


l-Cy. 


Because  of  the  regenerative  structure  of  A',  it  is  evident 
that  a  satisfies  the  equation  o  =  £[/li]  +  E[Ct]o.  Thus,  if 
E[Ct ]  <  1,  it  follows  that  a  is  finite  and  can  be  expressed 
as 


Hence,  the  infinite-horizon  discounted  cost  for  a  regenera¬ 
tive  process  can  be  expressed  in  terms  of  a  ratio  estimation 
problem;  see  Fox  and  Glynn  (1989)  for  further  details. 


EXAMPLE  3.  Let  X  be  a  regenerative  process  as  in 
Example  2,  and  assume  that  A'  has  right-continuous  paths 
with  left  limits.  Let  F  be  a  non-empty  subset  of  the  state 
space  S,  and  let  r(F)  =  inf { I  >  0  |  X(t)  €  F]  be  the  first 
hitting  time  of  the  subset  F.  Then, 


0  =  E[r(F)} 


is  the  mean  hitting  time  of  F.  Such  expectations  are  of 
interest,  for  example,  in  the  reliability  setting,  in  which 


case  r(F)  would  typically  correspond  to  the  system  failure 
time,  and  T(l)  to  a  time  at  which  the  system  is  brought 
back  to  an  “as  good  as  new”  state.  Let 

A,  =  min[r(£),  T(l)] 

By  =  7[r(F)  <  T(l)], 

where  /  denotes  the  indicator  function.  If  P[r(F)  <  oo]  > 
0  (note  that  this  is  equivalent  to  requiring  that  P[r(E)  < 
T(l)]  >  0),  it  is  easily  shown  that 


See  Goyal  et  al.  (1991)  for  additional  details.  Thus,  the 
mean  hitting  time  of  a  regenerative  process  can  be  formu¬ 
lated  in  terms  of  the  ratio  estimation  problem. 

EXAMPLE  4.  Let  X  be  a.  real-valued  random  variable 
and  let  C  be  an  event  with  P(C)  >  0.  Suppose  that  we 
wish  to  estimate 

<*  =  E[X  |  C], 

namely  the  conditional  expectation  of  X ,  given  that  the 
event  C  has  occured.  If  £[1X1]  <.  oo,  then  we  can  express 
or  in  terms  of  the  ratio  a  =  £[Ai]/£[£i],  where 

Ay  =  X  /(C) 

By  =  /(C). 


Hence,  conditional  expectations  are  expressable  in  terms 
of  the  ratio  estimation  problem. 

Thus,  the  ratio  estimation  problem  arises  in  a  variety 
of  different  applications  contexts.  We  shall  now  introduce 
a  decision  parameter  vector  9  into  the  discussion.  For  each 
9  €  HP*,  let  P»  be  the  probability  measure  associated  with 
the  parameter  value  9,  and  let  E»  be  its  corresponding 
expectation  operator.  In  addition,  we  shall  permit  the 
random  variables  A(0)  and  B(9)  to  depend  explicitly  on 
9  €  IR-.  Then,  for  each  9  6  IR*,  the  ratio  of  expectations 
can  be  expressed  in  the  form 


a(9)  = 


where  u(9)  =  Et[A(9)]  and  t(9)  =  Ea[B(9)].  Given  our 
above  examples,  computing  the  gradient  of  such  a  ratio 
o  (9)  is  useful  for  sensitivity  analysis  or  optimization  of 
any  of  the  following  :  steady-state  costs  or  rewards  in 
regenerative  processes;  infinite-horizon  discounted  costs; 
mean  time  to  failure  in  reliability  systems;  conditional  ex¬ 
pectations  and  probabilities. 


3.  CONFIDENCE  INTERVALS  FOR  GRADI¬ 
ENT  ESTIMATORS  OF  RATIOS 

Let  0O  €  IR'  be  fixed.  In  order  for  the  gradient  esti¬ 
mation  problem  to  make  sense,  we  shall  require  that  both 
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u()  and  ((■)  have  gradients  at  9  —  90.  We  shall  fur¬ 
ther  assume  that  there  exists  unbiased  estimators  for  not 
only  t»(0o)  and  t(9o),  but  also  their  gradients  Vii(0o)  and 
Vf(0o).  Focussing  now  on  the  i-th  component  of  the  gra¬ 
dient,  we  shall  specifically  assume  that  there  exist  jointly 
distributed  random  variables  (A,  B,C,D)  such  that 


E[A]  = 

v(90) 

E[B ]  = 

((So) 

E[C]  = 

d.u(fi0)‘= 

d 

d9, 

„(9) 

£[f>]  = 

d.((90)= 

d 

99. 

<(®)| 

where  d,  denotes  the  partial  derivative  with  respect  to  9,, 
and  9,  is  the  i-th  component  of  9. 

There  is  now  a  great  deal  of  literature  on  various 
ways  of  constructing  unbiased  estimators  for  d, u(9o)  and 
d,t(6o).  The  two  principal  approaches  that  have  been  ex¬ 
plored  are  likelihood  ratio  gradient  estimation  (see  Glynn 
1990  for  a  survey)  and  infinitesimal  perturbation  analysis 
(see  Glasserman  1991).  For  links  between  the  two  meth¬ 
ods  and  for  a  general  survey,  see  L'Ecuyer  (1990,  1991). 

We  shall  now  assume  that  it  is  possible  for  the  sim¬ 
ulator  to  generate  a  sequence  {(A,,  B,,  C,,  D,),  j  >  1} 
of  i.i.d.  replicates  of  the  random  vector  (A,  B,C,  D).  In 
each  of  the  problem  settings  described  in  Section  2,  this 
is  typically  straightforward . 

To  estimate 


d,a(90) 


e(60)d,v(90)  -  u(60)d,t(e0) 
F(9o) 

dm(9o)  —  o(90)d,t{9o) 

m  ’ 


the  natural  estimator  to  use  is 


where 


;=  i 

n 


;  =  > 


Our  first  proposition  states  that  under  reasonable  con¬ 
ditions,  4,(n)  is  a  consistent  estimator  for  d,o(0o).  The 
proof  is  straightforward  and  therefore  omitted. 

PROPOSITION  1.  Suppose  that  £[|.Ai  |  -f  |Bi |  +  |Ci|  -(- 
|Z>i|]  <  oo  and  that  E[B\]  0.  Then, 

lim  i,(n)  *=  3,o(0o)-  I 

n— .oo 

To  develop  a  confidence  interval  methodology  ford,(n), 
we  need  a  central-limit  theorem  (CLT)  for  the  estimator. 
Let 


Z,  =  A]  -  a(9o)B, 

W,  —  Cj  —  a(9o)D}  —  d,a(9o)Bj 

and  note  that  under  the  assumptions  of  Proposition  1, 
E[Zj ]  =  E[Wj ]  =  0.  This  observation  is  an  important 
element  in  the  proof  of  the  following  theorem. 

THEOREM  1.  Assume  that  E[Z\  +  W,2]  <  oo.  If,  in 
addition,  the  conditions  of  Proposition  1  are  in  force,  then 

\/n[Mn)  —  3,o(90)]  =>  crN (0, 1) 
as  n  — »  oo,  where 

,2  _  E[W,  -  (g[gi]/E[g1])Z^  _ 

(£[*«])’ 

Theorem  1  has  been  previously  established,  using  dif¬ 
ferent  methods,  by  Reiman  and  Weiss  (1989)  in  the  con¬ 
text  of  likelihood  ratio  gradient  estimation  for  regenera¬ 
tive  steady-state  simulation.  Their  expression  for  the  vari¬ 
ance  constant  <r2  is  formally  different,  but  algebraically 
identical. 

The  final  step  needed  to  develop  a  confidence  interval 
methodology  for  6,(n)  is  the  construction  of  an  appropri¬ 
ate  estimator  for  o1 .  Let 


v(n)  = 

where 

^-(Dn/Bn)2j}7 

j=»i 

(Bn)3 

=  Aj  —  otnBj 

=  Cj  —  otnDj  —  . 

The  next  proposition  gives  conditions  under  which  v(n) 
is  strongly  consistent  for  o3 .  The  proof  is  straightforward 
and  therefore  omitted. 

PROPOSITION  2.  Suppose  that  E[A \  +  Bf  +  C?  + 
D*]  <  oo.  If  E[Bi )  /  0,  then 


lim  r(n)  *=  <rJ.  I 


and 


a 


—  An/  Dn. 
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We  note  that  if  v(n)  is  computed  via  a  two-pass  ap¬ 
proach  in  which  ar„  and  S ,(n)  are  computed  in  the  first 
pass  through  the  data  { (-4^ ,  B,  .Cj ,  D}),  I  <  j  <  n]  and 
the  sum  of  squares  computed  the  second  pass,  then  it 
is  essentially  guaranteed  that  •)  will  be  computed  as 
a  non-negative  quantity  on  ant  nite-precision  computer. 
More  importantly,  this  means  of  computing  u(n)  is  likely 
to  be  more  stable  numerically  than  that  associated  with 
the  computation  described  in  Reiman  and  Weiss  (1989). 

We  are  now  ready  to  describe  a  general  confidence  in¬ 
terval  methodology  for  estimating  partial  derivatives  of 
ratios.  Suppose  that  we  wish  to  compute  a  100(1  -  6)% 
confidence  interval  (or  d,a($a).  We  use  the  following  pro¬ 
cedure  : 


Define 

W,(i)  =  Cj(i)  ~  «(*o )Dj(i)  -  d,a{80)B}. 

We  are  now  ready  to  state  a  joint  CLT  for  (a„,  5i(n),  . . . , 

M»))- 

THEOREM  2.  Assume  that  )+Z>?(l)4- 

•  ■  •  +  C?(d)  +  D\(d)]  <  oo.  //  £[Bi]  *  0,  then 

y/n[a„  -  or{0o),  Sj(n)~  d\  a[80), .... 

«d(f»)  -  dda(do)]  £[Bi]  =>  N(0,C) 

oj  n  -*  oo,  where  C  =  ( C,} ,  0  <  i,j  <  d)  is  a  covariance 
matrix  whose  elements  are  given  by 


Algorithm  CL 

1.  Generate  {(A},  B,,C},  Dj),  j  >  1). 

2.  Compute  a„  and  S,(n). 

3.  Compute  v(n)  (using  the  two-pass  approach  de¬ 
scribed  above). 

4.  Find  z(6)  such  that  P[N( 0, 1 )  <  z(6)]  =  1  —  5/2. 

5.  Compute 

Ln  =  5,(n)  -  z[S)\/v(n)/n) 

Un  =  5,(n)  +  z(b),Jv(n)/n).  | 

Then,  [£„,  £/„]  is  an  (approximate)  100(1—5)%  confidence 
Interval  for  d,a(0o)-  In  particular,  if  the  conditions  of 
Proposition  2  are  in  force  and  o2  >  0,  then 

Um  P[dia(90)  g  [L„,Un]]  =1-5. 

n-»  oo 

We  conclude  this  section  with  a  brief  discussion  of  the 
problem  of  generating  a  confidence  region  for  the  vector 
(a(6o),dia(8o), . . .  ,ddo(fi  )).  A  joint  confidence  region 
could  be  of  potential  interest  in  a  number  of  optimization 
settings,  since  virtually  all  iterative  (deterministic)  opti¬ 
mization  algorithms  choose  their  search  direction,  at  each 
iteration,  by  considering  the  full  gradient. 

Let  C(i)  and  D(i)  be  unbiased  estimators  for  d,u(8o) 
and  dil(8a),  so  that 

E[C(.)]  =  d,u(80) 

E[D(i))  =  d,f(9a). 

If  {(Aj.Bj.CjW.DAl) . C,(rf),D;(d)),  1  <  j  <  n} 

is  a  set  of  n  i.i.d.  replicates  of  the  random  vecto 
(A,B,C(  1),  29(1), . .  •  ,C{d),  29(d)),  then  the  estimators  o, 
5i(»),  5d(n)  can  be  constructed  from  the  sample  in 

the  obvious  way,  namely 

Of  n  — —  An/ Br, 

5i(n)  =  (Cn{t)-anDn(i))/Bn. 


Coo  =  E[2?) 

Co.  =  Cio 

a'  -  *[(-•«-»■) 

for  1  <  »,  j  <  d.  ■ 

The  proof  of  this  theorem  mirrors  that  of  Theorem  1 
and  is  therefore  omitted. 

A  procedure  for  producing  asymptotically  valid  con¬ 
fidence  regions  for  (a  (So),  dia(9o),  ...dda(8o))  can  now 
easily  be  derived,  using  the  same  ideas  as  those  described 
earlier  in  this  section  for  9,o(9o)-  We  leave  the  details  to 
the  reader. 


4.  LOW  BIAS  ESTIMATION  FOR  THE  GRA¬ 
DIENT  OF  A  RATIO 


Since  the  gradient  of  the  ratio  is  a  nonlinear  function 
of  the  expectations  E[A],  E[B],  E[C  'I,  2?[29(1)],  ..., 
E[C(d) ],  £[  29(d)],  it  follows  that  the  e  .ator  5<(n)  is,  in 
general,  biased  for  d,a(9o). 

We  will  now  proceed  to  (formally)  derive  a  bias  ex¬ 
pansion  for  {<(»).  The  proof  of  Theorem  1  shows  that 


6.(n)-di0(9o)  = 


Wn  -  ( Dn/Bn)Zn 

Bn 


(1) 


We  would  like  to  approximate  the  expectation  of  that.  We 
note  that  since  Bn  is  dose  to  p  =  B[Bt]  for  large  »,  we 
can  use  the  power  series  expansion  for  /(*)  =  (1  -  x)-1 
to  obtain 
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« 


2/1  —  Bn 

s  ' 

Using  this  approximation  in  (]),  we  find  that 
£.(n)  —  d,a(9o) 


2/i  —  B, 
n2 

l  1  2n  Dn  \ 

i,V"-  n.  ) 

2Wn  ( 

ZnDn  -  B„\V„ 

2  2nD„ 

M 

n2 

/‘Bn 

~  2lV"  , 

ZnDn  -  B„\Vn 

2ZnDn(2,L-  Bn) 

/» 

n2 

n3 

2Wn  _ 

ZnDn  +  Bn\Vn 

2  2n£)r\Bn 

whi-re  Bn  =  B„  —  /i.  Recall  that  E[Wj\  =  £[2j]  —  0. 
Observe  (hat  for  .  ft  j,  E[B,W,]  =  E{Bi]E[W}]  =  0,  since 
B,  and  W,  are  independent.  Therefore, 

E[B„Wn]  =  E[D,W})  = 

•  =1  ;=  1 

Similarly,  E[Z„D„]  =  E[ZlDl )/n.  Also,  E[Z,D}(Bk  - 
p)]  =  0  whenever  i  it  h.  Therefore, 

E[Z„D„B„]  =  g[g?P«(fli  -/»)] 

n* 

,  ("  -1  )E[Z,  (ft  -/Q]E[Q. 1 

=  ZMME1  +  0(1/B). 

Now,  taking  the  expectation  in  (2)  yields 
£[«.(*)]  -  d.o(0o) 

^  2E[£,Bi)E[D,)  -  ,tE[BlW1+ZlDl) 

TJ/l3 

This  bias  approximation  suggests  an  obvious  means  of  re¬ 
ducing  the  bias  of  gradient  estimators  for  ratios.  The  idea 
is  to  estimate  the  bias  term  and  correct  for  it  by  subtract¬ 
ing  off  the  estimated  bias.  In  this  case,  this  approach  leads 
to  the  estimator 


=  «.(»)  + 


2, d}) 

n'Dl 


where  Z}  and  fVj  are  defined  just  before  the  statement  of 
Proposition  2  in  Section  3. 


Under  appropriate  regularity  hypotheses,  and  by  ap¬ 
plying  techniques  similar  to  those  used  in  Glynn  and  Hei- 
delberger  (1991),  one  can  rigorously  prove  that  i,(n)  re¬ 
duces  the  asymptotic  bias,  in  the  sense  that 

£[6,(n)]  =  d,a(0o)  +  o(l/n). 

A  second  approach  that  is  frequently  used  to  correct 
for  “nonlinearity  bias’  of  the  above  type  is  to  “jackknife” 
the  estimator.  Specifically,  for  1  <  j  <  n,  let 

n 

E  ■<* 

*=1. 

“"())  =  - n - 

E  * 

*=1,  k^tj 

(  it,  Cr)  -  a«0)  it  D> 

a  \*=1.  /  fcssl,  k^lj 

pn{j)  =  - - - 

E 

*«<»  =  ”a„ (»  -  (n  -  1  )0n(j). 

Then, 

6'  (n)  = 

j=i 

is  the  jackknife  estimator  for  d,a(90).  Also, 

3  Kn) 

where 

(,J(n))J  =  TrTD^)-^(n))J 

:=i 

is  a  consistent  variance  estimator.  As  in  the  case  of  the 
estimator  l.(n),  one  can  prove  rigorously  (under  suitable 
regularity  hypotheses)  that  the  estimator  S/(n)  reduces 
asymptotic  bias,  in  the  sense  that 

£[«/(«)]  =  dia(90)+o(l/n). 

It  turns  out  that  the  improved  bias  characteristics  of 
these  estimators  are  costless  relative  to  the  variance,  in  the 
sense  that  the  estimators  l,(n)  and  6f(n)  obey  precisely 
the  same  CLT  as  does  6<(n).  Hence,  the  estimators  exhibit 
the  same  degree  of  asymptotic  variability. 

THEOREM  3.  Assume  that  E[A3  +  B2  +  C2  +  D2]  <  oo 
and  that  £[Bi)  ft  0.  Then, 

Vn(Si(n)  -  d,a(8o))  =*  «rAf(0, 1) 
y/n(6f{n)  -  d,a(0o))  =*•  <rN(0, 1) 

where  a3  is  the  same  constant  as  in  Theorem  1.  | 


Table  1  Bias,  Half-Widths,  and  Coverages  of  Confidence  Intervals  for  8  =  0.2 
_ (o(tf)  =  0.25  and  a‘(9)  =  1.5625) _ 


n  =  10 

n  =  100 

n  =  1000 

bias 

half-width 

cover. 

bias 

half-width 

cover. 

bias 

half-width 

cover. 

-.010±.002 

.129±.002 

,77±.01 

-.001±.001 

.059±.001 

.90±.01 

-.OOOi.OOl 

.021±.001 

94±.01 

fm 

-,001±.002 

,149±.002 

,80±.01 

-,000±.001 

.061±.001 

.90±.01 

-OOOi.OOl 

.021±.001 

.94±.01 

-.401±.024 

,914±.021 

.43±.01 

-113±.013 

.857±,014 

.74±.01 

-004±.005 

,391±.004 

,88±.01 

£ m 

-,120±.035 

1.345±.038 

.52±.01 

-,014±.013 

.921±.017 

.76±.01 

-001±.005 

.395±.004 

,88±.01 

Table  2:  Bias,  Half- Widths,  and  Coverages  of  Confidence  Intervals  for  8  =  0.5 
_  (a(6)  =  1.0  and  o’(fl)  =  4.0) _ 


n  =  10 

n  =  100 

n  =  1000 

bias 

half-width 

cover. 

bias 

half-width 

cover. 

bias 

half- width 

cover. 

f  ‘  ‘ 

On 

-.135±.009 

,440±,005 

,57±,01 

-,019±.004 

,313±.004 

.82±.01 

-  003±.002 

.125±.001 

,92±.0l 

-036±.011 

.534±.0 12 

.65±.01 

-001±.004 

,339±.004 

.84±.01 

-001±.002 

.126±.001 

.92±.01 

S(n) 

-2.072±.050 

1.692±.057 

,26±.01 

-,407±.048 

2.333±.051 

■57±.01 

-056±.019 

1.493±.022 

.81±.01 

8J(n) 

-1.152±.089 

2.694±.084 

,37±.01 

-.049±.063 

2.919±.077 

.63±.01 

-012±.020 

1.541±.022 

.82±.0l 

5.  A  NUMERICAL  EXAMPLE 

We  will  now  illustrate  some  of  the  ideas  developed  in 
this  paper  with  a  numerical  example.  We  consider  the 
steady-state  : ojourn  time  of  a  customer  in  an  M/M/1 
queue  with  arrival  rate  A  =  1  and  mean  service  time  8, 
where  0  <  8  <  1.  The  sojourn  time  is  the  sum  of  the 
time  spent  by  a  customer  waiting  for  service,  plus  that 
customer’s  service  time.  Let  A'o  =  0  and,  for  »  >  1,  let  X, 
be  the  sojourn  time  of  the  i-th  customer,  starting  from  an 
empty  system.  It  is  well  known  that 

Xi  .=  (X._,  -r/.)+  +0C., 

where  {iq  ■  •  -  ■}  is  a  sequence  of  i.i.d.  exp(l)  ran¬ 

dom  variables.  For  this  model,  the  mean  steady-state  so¬ 
journ  time  or(0)  can  be  computed  in  closed  form: 

a(8)  =  9/(l-8). 

Hence, 

<*'(*)  =  1/(1  -  of 

This  system  regenerates  when  customers  arrive  to  an 
empty  queue.  Consequently,  as  discussed  in  Example  1  of 
Section  2,  the  steady-stale  mean  sojourn  lime  of  a  cus¬ 
tomer  can  be  expressed  in  terms  of  a  ratio  estimation 
problem,  and  the  methodology  of  this  paper  is  therefore 
applicable.  It  is  also  straightforward  to  apply  the  likeli¬ 
hood  ratio  method  for  gradient  estimation  to  this  problem 


(see  L’Ecuyer,  Giroux,  and  Glynn  1991),  thereby  obtain¬ 
ing  the  required  unbiased  estimators  for  the  numerator 
and  denominator  of  the  ratio  (as  discussed  in  Section  3). 
It  turns  out  that  while  infinitesimal  perturbation  analysis 
can  be  applied  to  obtain  strongly  consistent  steady-state 
gradient  estimators  for  this  problem,  it  fails  to  give  unbi¬ 
ased  estimators  of  the  gradient  of  the  numerator  and  of  the 
gradient  of  the  denominator  of  the  regenerative  ratio  for¬ 
mula.  See  Heidelberger  et  al.  (1988)  for  further  details,  as 
a  consequence,  the  theory  of  this  paper  is  not  applicable  to 
the  infinitesimal  perturbation  analysis  steady-state  gradi¬ 
ent  estimator  for  this  problem.  But  on  the  other  hand,  the 
infinitesimal  perturbation  analysis  derivative  estimator  is 
itself  the  estimator  of  a  ratio  of  expectations,  so  that  one 
can  apply  the  standard  theory  relative  to  the  construction 
of  confidence  intervals  for  ratios  of  expectations  (Iglehart 
1975,  Wolff  1989). 

Tables  1  and  2  report  the  experimental  results  ob¬ 
tained  for  this  example.  Simulation  runs  were  carried  out 
at  two  parameter  values,  namely  0  =  0.2  and  $  =  0.5, 
using  n  =  10,  100,  and  1000  regenerative  cycles.  A  to¬ 
tal  of  four  estimators  were  considered  in  this  experiment, 
namely  the  ratio  estimator  arn  for  a(0)  and  its  jackknifed 
analog  ai  (see  Iglehart  1975),  and  the  derivative  estima¬ 
tor  5(«)  (Section  3)  for  <*'(0)  and  its  corresponding  jack¬ 
knifed  analog  8J(n)  (Section  4).  Standard  regenerative 
confidence  intervals  were  constructed  for  the  estimator 
On,  and  the  confidence  interval  methodology  of  Section 
3  was  used  to  analyze  0(n).  For  the  jackknifed  versions, 
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confidence  intervals  were  constructed  based  on  the  vari¬ 
ance  estimator  (sJ(n))2  given  in  Section  4.  At  each  of  the 
two  parameter  values  and  three  choices  of  n  (the  number 
of  regenerative  cycles),  a  total  of  10,000  95%  confidence 
intervals  was  replicated  for  each  of  the  four  estimators. 
From  that,  we  are  able  to  report  estimates  for  the  bias,  ex¬ 
pected  half-width,  and  coverage  (the  probability  that  the 
quantity  being  estimated  lies  in  the  confidence  interval), 
as  well  as  95%  confidence  intervals  for  the  bias,  expected 
half-width,  and  coverages  themselves. 

One  can  observe  that  for  small  n,  for  all  estimators, 
the  coverage  is  really  lower  that  what  is  to  be  expected. 
This  bad  behavior  gets  worse  when  9  increases  (heavier 
traffic).  Jackknifing  clearly  reduces  the  bias  significantly. 
It  also  gives  a  better  coverage  for  small  n,  but  usually  at 
the  expense  of  a  wider  confidence  interval.  For  small  n, 
the  coverage  is  too  low  anyway.  For  larger  «,  jackknifing 
still  helps  reducing  the  bias,  but  (perhaps  surprisingly) 
does  not  improve  the  coverage  significantly.  Of  course, 
this  is  just  a  particular  illustration,  and  one  must  be  care¬ 
ful  about  drawing  any  general  conclusions  from  these  nu¬ 
merical  results. 

6.  CONCLUSION 

Ratio  estimation  problems  arise  in  many  different  ap¬ 
plications  settings.  When  estimation  is  to  be  used  to  an¬ 
alyze  the  sensitivity  of  (or  to  optimize)  a  system  in  which 
the  ratio  estimation  problem  occurs,  the  results  of  this 
paper  become  pertinent.  We  have  derived  a  numerically 
stable  confidence  interval  procedure  for  computing  partial 
derivatives  of  such  ratios,  and  have  developed  the  appro¬ 
priate  joint  CLT’s  necessary  to  extend  this  methodology 
to  the  computation  of  confidence  regions  for  the  full  gradi¬ 
ent  of  the  ratio.  In  addition,  we  have  discussed  low-bias  es¬ 
timators  for  computing  such  partial  derivatives.  We  have 
also  described  preliminary  computational  experience  with 
some  of  the  methods  developed  in  this  paper. 
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APPENDIX 

Proof  of  Theorem  1.  We  note  that 

5„[«,(n)  -fto(tfo)] 

=  Cn  ~  OnDn  -  d,o{90)Dn 
-  W„  -  (o„  -  n(90))D„ 

=  Wn  -  (D„/Bn)Z„ 


=  Wn  -  (E{D\)f E[B\))2n 

-  ( Dn/Bn  -  £[£>,]/£[£,])  Zn. 

Clearly,  *Jn2n  =>  (£(Z?])>/2 N(0, 1)  as  n  —  oo  and 
DnlB„  —  E[D\]/ E[B\]  as  n  —  oo.  It  follows,  by  the 
converging-together  principle,  that 

>/£  (Dn/Bn  -  £(D,  ]/£[£,])  Zn  =>  o 

as  n  — >  oo.  The  CLT  for  i.i.d.  random  variables  also 
proves  that 

sft(Wn-  (E[D,]/E[B,])Zn)  =>  E[Bi]*N(0, 1) 

as  n  — -  oo.  A  second  application  of  the  converging- 
together  principle  then  yields 

VZB„  («,(»)  -  3,a(90))  =>  E[Bt]<rN( 0, 1). 

One  final  application  of  the  converging-together  principle 
(note  that  B„  E[B i]  as  n  —  oo)  proves  the  theorem. 

■ 
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